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Grid-Free  Electromagnetic  Particle  Simulations 

Andrew  J.  Christlieb 
Department  of  Mathematics 
Michigan  State  University 

Abstract  Basic  plasma  science  plays  an  increasingly  significant  role  in  applications  of  impor¬ 
tance  to  the  United  States  Air  Force.  Many  of  these  applications  require  a  fully  kinetic  description 
in  at  least  part  of  the  domain.  The  most  common  approach  is  to  use  a  fully  Lagrangian  framework, 
where  the  model  is  reduced  to  tracking  the  evolution  of  test  particles  in  phase  space.  Of  the  many 
varieties,  the  most  accepted  approach  is  Particle-In-Cell  [1].  The  PI,  with  his  collaborators,  is  de¬ 
veloping  an  alternative  approach,  the  Boundary  Integral  Treecode  (BIT)  [2].  BIT  is  based  on  fast 
summation  algorithms  and  boasts  arbitrary  accuracy.  In  a  range  of  numerical  experiments,  the  all 
scale  resolution  of  BIT  has  proven  to  provide  a  substantial  improvement  over  PIC  [2].  A  major 
objective  in  this  currant  work  is  to  develop  a  grid-free  electromagnetic  formulation  of  BIT. 

Summary  Over  the  past  four  years,  while  developing  our  methodology  for  a  grid-free  electro¬ 
magnetic  formulation  of  BIT,  several  significant  topics,  which  need  to  be  addressed,  have  become 
central  in  our  research  efforts.  These  include  the:  1)  Extended  our  ID  analysis  of  the  impact  of  ker¬ 
nel  regularization  on  the  accuracy  of  the  solution  in  the  vicinity  of  boundaries  to  3D  (in  preparation 
[8]-  Dr.  Christlieb,  Dr.  Cartwright  and  Dr.  Ong),  2)  worked  on  development  of  Boundary  Integral 
Corrected  PIC  (in  preparation  -  Dr.  Christlieb,  Dr.  Cartwright  and  Dr.  Ong),  3)  a  multi-scale  ap¬ 
proach  to  time  stepping  of  coupled  stiff  systems  based  on  our  high  order  semi-implicit  integrators 
[3,4,5],  4)  applied  our  parallel  time  integrator  [6,7]  to  problems  that  to  time  stepping  of  one  million 
particles  using  3D  BIT  to  compute  fields  and  demonstrated  that  even  with  the  problem  does  not  fit 
in  cash  we  still  get  4th  order  in  wall  clock  time  of  forward  Euler  (in  preparation  -  Dr.  Christlieb, 
Dr.  MacDonald  and  Dr.  Ong),  5)  development  of  arbitrary  order  asymptotic  preserving  methods 
for  multi-scale  problems  (under  development  -  Dr.  Christlieb),  6)  development  of  arbitrary  order 
adaptive  methods  based  on  explicit  and  semi-implicit  integral  differed  correction  (work  complete  - 
in  preparation  -  Dr.  Christlieb,  Dr.  MacDonald,  Dr.  Ong),  7)  a  study  of  the  impact  of  regularization 
on  solutions  to  the  Schrodinger  equation  -  (in  preparation  -  Dr.  Christlieb,  Dr.  Krasny  and  Dr.  Ong), 
8)  translation  and  adaptation  of  a  3D  c-i-i-  BIT  code  to  parallel  CUDA  for  GPU  architecture  (under 
development  -  Dr.  Christlieb,  Dr.  Johnston  and  Dr.  Ong),  9)  development  of  implicit  time  step¬ 
ping  methods  based  for  Maxwell’s  equations  using  the  BIT  metrology  (currently  modifying  3D 
BIT  code  to  handle  3D  Maxwell  kernel  -  Dr.  Christlieb,  Mr.  VanGroningen  and  Dr.  Ong),  10) 
extension  of  BIT  to  radiative  transport  problems  ([9]/method  under  development  -  Dr.  Christlieb, 
Dr.  Hitchon,  Dr.  Eawler  and  Dr.  Eister)  11)  extension  to  2D  of  an  adaptive  mesh  refinement  frame¬ 
work  based  on  WENO  (Dr.  Christlieb,  Dr.  Qiu  and  Dr.  Shen  -  Goal  is  AMR  VP  solver  based  on 
WENO)  [10].  All  of  these  topics  are  linked  by  a  desire  to  bridge  the  multiple  spatial  and  temporal 
scales  which  arise  in  plasma  problems.  We  now  give  a  brief  summer  of  the  progress  on  our  time 
stepping  methods,  analysis  of  BIT,  extension  of  BIT  to  Maxwell’s  equations  and  our  high  order 
WENO  AMR  framework.  At  the  end  of  this  document,  I  have  attached  our  relevant  publications. 
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Integral  Deferred  Correction  Methods  The  Integral  Deferred  Correction  (IDC)  framework 
[3,4]  being  developed  by  the  PI  and  his  collaborators,  has  been  shown  as  one  possible  way  to 
systematically  generate  high  order  RK  methods  based  on  a  low  order  RK  method.  IDC  belongs 
to  the  family  of  integrators  known  as  defect  correction  algorithms.  Given  a  system  of  the  form 
^  =  /(t,y)  with  y(to)  =  yo,  IDC  works  as  follows:  let  v®  be  an  approximation  to  y(t)  given  by  a 
p^*  order  RK  method  on  the  time  steps  to,  ■  •  • ,  let  v(t)  be  a  polynomial  of  degree  at  most  m  pass¬ 
ing  thought  the  V®;  define  the  error  ase{t)  =y{t)  —v{t)  and  the  residual  as  r(t)  =  —  f{t,v{t)); 

let  Q{t)  =  e{t)  +  r{x)dx;  taking  the  derivative  of  the  equation  for  e{t)  and  re-writhing  the  equa¬ 
tion  in  terms  of  2 (t),  we  arrive  at  the  differential  form  of  our  IDC  corrector  equation, 

^  =  /(^^e  +  Vo  +  ^  /(x,v(x))Jx^  -/(t,v(t))  . 

The  corrector  equation  is  solved  to  correct  v®  differs  from  other  defect  correction  methods  in  that 
the  residual  does  note  appear  in  differential  form,  rather  the  integral  of  the  residual  is  built  into  the 
^  equation.  This  inherently  makes  the  method  more  robust  than  standard  defect  correction.  ^ 
is  solved  with  the  same  order  RK  method  and  used  to  compute  v|  <-v^  +  Qi  -  fi‘  r{x)dx.  The 
PI,  along  with  his  current  and  former  post  docs  Dr.  Ong  and  Dr.  Qiu,  have  show  that  under  cretin 
assumptions,  r  —  1  successive  approximations  of  the  corrector  equation  solved  with  a  order 
explicit  RK  method,  results  in  a  r  x  order  method,  up  to  a  maximum  order  L,  which  is  only 
determined  by  the  accuracy  in  the  approximation  of  the  integral  residual,  r{x)dx.  In  practice, 
we  have  used  a  p^^  order  implicit  and  semi-implicit  RK  in  the  predictor  and  corrector  and  observed 
the  same  result.  Dr.  Morton,  together  with  the  PI,  Dr.  Ong  and  Dr.  Qiu,  has  extended  the  theoret¬ 
ical  results  in  [3,4]  to  both  implicit  and  semi-implicit  RK  methods.  The  PI  is  currently  working 
on  extending  these  ideas  to  class  of  semi-implicit  RK  methods  know  as  Asymptoticly  Preserving 
(AP)  methods.  AP  methods  are  a  systematic  philosophy  for  time  stepping  of  systems  with  stiff 
relaxation.  Given  a  system  with  a  small  parameter,  say  the  Knudson  number  in  the  collisional 
Boltzmann  equation,  AP  methods  are  a  two  step  proceeder  for  developing  numerical  methods  that 
give  uniform  convergence  to  the  limiting  system  as  the  small  parameter  goes  to  zero,  regardless  of 
the  time  step.  Step  one  involves  either  mathematical  modeling  or  a  clever  re-casting  of  the  system 
in  form  where  taking  the  limit  as  the  small  parameter  goes  to  zero  yields  the  ‘correct’  limiting 
behavior.  For  the  Boltzmann  equation  of  neutral  gas  dynamics,  this  involves  scaling  the  system 
so  that  the  diffusion  limit  is  obtained  as  Knudson  number  goes  to  zero.  Step  two  involves  the 
development  of  a  consistent  Semi-Implicit  Runge-Kutta  (SIRK)  schema  designed  to  give  uniform 
convergence  as  the  small  parameter  goes  to  zero.  Then  implicit  part  of  the  SIRK  method  handles 
the  stiff  term,  in  this  example  the  collision  term,  implicitly  and  the  non  stiff  term,  here  transport, 
explicitly.  The  PI  has  show  that  the  IDC  framework  can  support  the  type  of  SIRK  that  shows  up  in 
an  AP  formulation.  It  remains  to  be  shown  that  the  new  method  that  results  from  embedding  the 
AP  SIRK  in  an  IDC  framework  will  again  be  an  AP  method.  We  are  in  the  processes  of  applying 
the  IDC  AP  SIRK  to  a  range  of  test  problems  to  provide  numerical  evidence  that  the  resulting 
method  is  indeed  AP,  we  will  then  focus  on  rigorously  showing  this  to  be  the  case. 

We  have  also  focused  on  extending  our  method  to  a  parallel  time  integrator.  IDC  differs  from 
traditional  defect  correction  in  that  the  integral  from  of  the  residual  is  used  in  the  formulation  of 
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the  solution,  instead  of  the  differential  form  of  the  residual.  On  multi-core  architectures,  defect 
correction  methods  become  extremely  attractive  because  each  correction  step  can  be  decoupled 
from  the  prediction  and  prior  correction  steps.  In  so  doing,  the  methods  are  able  to  achieve  high 
order  accuracy  in  the  wall  clock  time  equivalent  to  that  of  the  prediction  step,  provided  multiple 
cores  are  used  for  the  computation.  We  refer  to  this  new  class  of  defect  correction  as  Revisionist 
IDC  (RIDC),  which  allow  for  the  corrections  to  be  computed  in  parallel.  The  parallel  integrator 
we  have  developed  is  the  first  integrator  to  leverage  multi-core  cpus  such  that  the  method 
achieves  order  in  the  wall  clock  time  equal  to  that  of  a  single  forward  Euler  step. 


BIT  Corrected  PIC  The  use  of  time  depended  Green’s  functions  in  a  fully  Lagrangian  frame¬ 
work  raises  the  issue  of  needing  to  track  a  time  history  for  each  test  particle.  One  proposed  ap¬ 
proach  to  controlling  the  need  for  a  time  history  is  to  make  use  of  fixed  course  mesh.  As  a  first 
step  towards  the  electromagnetic  case,  we  are  exploring  the  use  of  a  coarse  fixed  mesh  for  domain 
decomposition  in  electrostatic  problems,  where  BIT  is  used  as  a  sub-cell  method  within  each  PIC 
cell.  The  methods  presented  here  differ  from  Particle-Particle-Particle-Mesh  (P3M)  in  that  local 
boundary  integrals  are  used  within  each  cell  to  provide  an  accurate  description  of  the  local  fields 
within  a  mesh  cell  [6].  We  demonstrate  that,  in  ID,  regularized  BIT  corrected  PIC  substantially 
reduces  numerical  heating,  even  when  Kx^Xd-  Because  of  the  negative  impact  regularization  has 
on  the  solution  field  near  boundaries,  we  have  developed  localizations  of  the  regularization  based 
on  Taylor  expanding  the  ID  non-regularized  greens  function  bout  the  regularized  kernel  and  have 
further  proven  that  these  are  convergent  expansions.  This  provides  a  possible  adaptive  approach 
for  dealing  with  regularized  kernels  near  fixed  non-periodic  boundaries.  Further,  we  expect  that 
Ar  should  increase  the  efficiency  of  large  scale  PIC  calculations,  enhancing  the  capability 
of  legacy  PIC  codes  by  adding  both  fidelity  (less  numerical  heating)  and  efficiency  (allowing  more 
flexibility  in  meshing  at  a  given  level  of  fidelity).  For  a  wide  class  of  electromagnetic  problems, 
where  PIC  is  used,  this  electrostatic  work  is  beneficial,  since  for  stability,  the  mesh  spacing  is  suf¬ 
ficiently  small  that  the  system  looks  electrostatic.  In  collaboration  with  AFRL/RDHE,  the  PI  is 
working  on  incorporating  this  into  ICEPIC.  Below  we  review  the  analysis  of  ID  BIT.  The  analysis 
has  just  recently  been  extend  to  the  2D  and  3D  case. 


ID  BI-PIC  -  Regularization  Near  cell  boundaries,  the  solution  of  regularized  BIT  diverges  from 
the  true  solution.  To  understand  this  error,  we  Taylor  expand  the  non-regularized  free  space  Green’s 
function  about  the  regularized  one.  We  establish  that  this  expansion  is  rigorously  convergent  as 
the  number  of  terms  goes  to  infinity  and  used  this  result  to  formally  establish  error  bounds  at  a  cell 
boundary.  In  addition,  this  expansion  has  a  compact  form  in  ID  and  gives  a  way  of  localizing  the 
regularization.  The  convergent  Taylor  expansion  of  {x  —  yY  -\- cP-  —  cP-  about  {x  —  yY  + 

(P, 


1  /  ^ 

\i=0 


(-l)'+i(2/-3)!! 
¥i\ 


-py  [{x-yY  +  P) 


.2^  (1-29/2 


(1) 


where  (  — l)'+^(2z  —  3)!!  =  (1)(1)(— 1)(— 3)(— 5) . .  .(7  —  2z)(5  —  2z)(3  —  2z),  shows  up  as  part  of 
the  constant  in  the  P  derivative  in  the  above  Taylor  expansion.  Using  a  similar  procedure  for 


3 


Mathematics 


AJ.  Christlieb 


dx 


gives, 


Final  Report  on  AFOSR  Grant  Number  FA9550-07-01-0092 


dG%{x\y) 

dx 


¥i\ 


{-d^)\x-y)  {{x-y)^  +  d^) 


-(2!'+1)/2 


(2) 


where  ( —  1 )' (2/  —  1 ) ! !  =  ( 1 )  (— 1 )  ( — 3) . . .  (5  —  2/)  (3  —  2/)  ( 1  —  2/) .  The  error  bounds  for  the  integral 
equation  with  regularized  field  for  a  cell  with  endpoints  at  a  and  (3,  with  potential  fields  given  at 
these  points,  using  an  N  term  approximation  for  can  be  written  as. 


d  /^d 


(^) 


-(27V+l)/2 


dy+ 


!£iP{y)TxGwiAy)dy- 


(3) 


Y.Wip{xi)^^G%{x\xi 


where  ^  G  ((v  —  y)^,  —  +  In  equation  (3),  the  term  (^) 


-(2Af+l)/2' 


is  strictly  deceasing 


on  the  interval  ^  G  {{x  —  kY,{x  —  kY-\-d'^)  and  has  its  maximum  value  and  slope  at  ^  =  (x  —  y)^. 
Hence  equation  (3)  is  bounded  above  by. 


TPyy  ^  Pmax  ^(-^+2)  j2(N+1) 
-  40rr(iV+2)^ 

d  ^d 
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APv+T)72 
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!Q.p{y)Tx^wiAy)dy  -  Y^wip{xi)^^G\i^{x\xi) 


(4) 


where  pmax  is  the  maximum  density  in  the  domain  Q.  and  a  and  |3  are  the  end  points  of  the  domain 
Q..  Equation  (4)  agrees  with  our  the  numerical  results,  predicting  the  correct  upturn  at  the  left 
boundary  and  the  correct  down  turn  at  the  right  boundary.  Further,  equation  (4)  demonstrates  that 
the  more  localized  the  regularization,  the  closer  the  defect  will  be  localized  to  the  boundary,  i.e., 
large  N  will  localize  the  error.  Figure  1  a)  shows  using  six  terms  of  the  equations  (1)  and  (2).  It 
is  clear  that  the  additional  terms  cause  the  jump  in  PIC-BIT-5  to  dramatically  decrease.  In  addition, 
the  interior  of  the  cell  the  PIC-BIT-5  and  PIC-BIT-5-ext  are  in  very  close  agreement.  Further,  using 
6  terms  in  the  approximation  of  Gj^,  picks  up  more  of  the  features  of  the  exact  solution  given  by 
non-regularized  BIT.  Figure  I  b)  shows  the  exact  solution  (blue  diamonds)  and  PIC-BIT-5  with  one 
(magenta  triangles),  two  (green  circles)  and  three  (red  dots)  terms  in  the  approximation  of  G®^.  It 
is  important  to  note  that  the  more  localized  the  approximation  to  G®^,  the  smaller  the  time  step 

d(jd 

will  need  to  be  in  the  numerical  simulation  so  as  to  resolve  the  transition  in  Our  work  shows 

ax 

dCj^ 

that  the  better  an  approximation  to  the  sharper  the  transition  that  needs  to  be  resolved.  So 
the  tradeoff  will  be  between  efficiency  and  localization  of  G‘^ . 


Time  Implicit  BIT  Building  on  the  idea  of  using  the  ‘fixed’  mesh  to  support  a  time  history,  the 
PI  and  Dr.  Ong  have  developed  a  method  of  lines  transpose  approach  to  wave  equations.  In  this 
case,  the  mesh  is  not  a  standard  finite  difference  mesh,  but  rather  a  way  of  tracking  the  time  evo¬ 
lution  of  waves  not  in  the  proximity  of  Lagrangian  test  particles.  The  mesh  points  do  not  have  to 
have  any  uniformity  to  them,  and  could  be  thought  of  as  ghost  particles,  which  could  be  done  with 
a  moving  mesh. 
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a)  b) 

Figure  1:  Foree  Vs.  Loeation:  a)  PIC  is  the  blaek  dashed  line,  BIT  is  the  green  line  with  dots,  BIT 
eorreeted  PIC  with  d  =  0  is  the  open  red  eireles,  six  term  BIT  eorreeted  PIC  with  d  =  ISx/ IQ 
is  the  blue  diamonds  and  extended  eell  six  term  BIT  eorreeted  PIC  is  the  magenta  triangles, 
b)  PIC  is  the  blaek  dashed  line,  BIT  is  the  blue  diamonds,  one  term  BIT  eorreeted  PIC  with 
d  =  Ax/ 10  is  the  magenta  triangles,  two  term  BIT  eorreeted  PIC  with  J  =  Ax/ 10  is  the  green 
eireles,  and  three  term  G^  BIT  corrected  PIC  with  d  =  Ax/ 10  is  the  red  dots. 


Time  Implicit  BIT  for  Wave  Type  Equations:  The  key  idea  is  to  consider  the  transpose  of  the 
standard  method  of  lines  methodology,  i.e.,  we  choose  to  discretize  in  time  and  directly  solve  the 
resulting  Helmholtz  equation  using  an  integral  formulation.  Observe  that  Maxwell’s  equations  can 
be  cast  as 

-^Ett  —  y^E  ~  —lUoJt - V p  ,  ~  y^B  =  /ToV  X  J  , 

where  p  =  /  fdv  is  the  charge  density,  J  =  fvfdv  is  the  current  density,  (E^,  Po)  are  the  permittivity 
and  permeability  and  c  is  the  speed  of  light.  To  illustrate  the  implicit  treecode  time  stepping 
methodology  for  E  and  B,  it  suffices  to  consider  the  wave  equation,  Utt  —  k^y^u  =  0  .  Using  a 
centered  difference  approximation  to  Ua  and  evaluating  at  time  level  n+l  gives. 


2„«+l 


k^At^ 

The  integral  solution  for  is 


y^u 

,n+l 


1 


U"+^  = 


k^At^ 


(-2m"  +  m"-1)  . 


,«+i 


(x)  = 


Q. 


k^At^ 


Gix\y)dTly+i  (m"+^VG- GVm”+M -nJ^, 

J  do. 


where  G(x|y),  the  free  space  Green’s  function  for  the  Helmholtz  operator  l{-)  =  (V^  — 

is  G(x|y)  =  ]^3  f  —  ||;k-_  j||2  and  yis  the  normalization.  The  volumetric 

and  boundary  integrals  are  approximated  at  the  midpoints,  and  the  resulting  sums  can  be  computed 
using  fast  summation  algorithms.  For  the  wave  equation,  the  method  is  able  to  take  time  steps  much 
larger  than  the  imposed  CFL  restrictions  for  an  explicit  integrator.  We  have  applied  our  approach 
to  the  wave  equation  in  ID,  2D  and  3D  with  Direcltet  boundary  conditions.  Additionally,  this 
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implicit  time  stepping  methodology  has  been  extended  to  higher  order  approximations  in  time. 
Some  interesting  behaviors  can  be  observed  in  the  simulations  below.  In  (a),  the  first  order 
scheme  converges  slowly  to  the  analytic  solution  while  in  (b),  the  third  order  scheme  converges 
much  quicker,  albeit  with  more  oscillations.  A  convergence  study  is  shown  in  (c).  In  this  example 
the  number  of  intervals  used  in  the  mid-point  approximation  to  the  volumetric  integral  is  M  =  4000. 
The  initial  condition  is  set  to 

Uo{x)  =  exp(— 

The  time  step  used  in  the  the  implicit  updates  are  At  =  with  N  E  {1,2,.  ..,5}.  Here  the  final 
time  is  tf  =  5.  For  all  results,  the  analytic  solution  was  used  to  start  the  time  marching  method.  In 
each  of  the  plots,  blue  curve  is  the  exact  solution,  the  green  curves  are  numerical  solutions  with 
At  E  {^5  ^5  and  the  red  curves  are  numerical  solutions  with  At  =  ^.  In  all  the  plots  we  see 
convergence  to  the  true  solution  as  At  — )■  0.  Observe  that  in  all  cases  At  is  much  larger  than  a  cell 
size.  Note  that  while  in  the  second  and  third  order  examples  display  oscillations  for  large  At,  the 
oscillations  rapidly  deceases  as  At  approaches  zero,  oscillations  are  not  visible  for  At  =  i. 


1st 

2nd 

a 

-  3rd 

* 

• 

(a)  first  order  time  differencing  (b)  third  order  time  differencing 


(c)  convergence  as  At  — 0 


Fast  Summation:  To  overcome  the  0{N^)  issue,  BIT  makes  use  of  a  fast  summation  algorithm. 
There  are  many  variants  of  fast  summation  algorithms,  in  this  case  we  are  making  use  of  the 
Treecode  algorithm.  The  idea  is  to  approximate  the  long  range  interaction  of  each  distant  cluster 
of  particles  as  a  moment  expansion  about  the  center  of  each  cluster.  Particle  clusters  are  creates 
with  a  hierarchal  tree  sorting  of  the  particles.  Fields  are  computed  using  a  recursive  divide  and  con¬ 
quer  approach,  which  makes  use  of  the  tree  [2] .The  method  is  0{NlogN),  a  substantial  speedup 
over  direct  summation.  (Note  that  PIC  is  0{M\ogM)  where  M  is  the  number  of  mesh  elements.) 

In  the  work  by  Lee,  Johnston,  and  Krasny  (JCP  09),  a  recursion  relation  is  developed  for  kernels 
of  the  from  Yexp(— cr)/r.  This  kernel  is  precisely  the  from  of  the  free  space  Green’s  function 
that  arises  in  the  3D  time  implicit  formulation  Maxwell’s  equations.  The  recursion  relation  is  a 
simple  modification  of  the  recursion  relation  used  in  electrostatic  BIT  and  are  in  the  processes  of 
modifying  our  3D  c-i-i-  electrostatic  BIT  code  so  that  we  will  be  able  to  quickly  develop  a  3D  time 
implicit  Maxwell  solver.  An  essential  step  is  to  consider  the  development  of  non-oscillatory  time 
stepping  methods  which  are  higher  than  first  order  in  time  for  the  time  implicit  Maxwell  solver. 

High  order  AMR  based  on  finite  difference  WENO  As  a  side  project,  we  have  been  working 
seedily  to  develop  a  high  order  strategy  for  Adaptive  Meth  Refinement  based  on  Weighted  Essen- 
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daily  Non-Oscillatory  finite  difference  methods.  The  new  method  is  the  first  which  is  3rd  order  in 
time  and  5th  order  in  space.  The  key  ideas  that  have  been  introduced  which  differ  from  standard 
AMR  are  that  the  method:  1)  treat  refinement  patches  that  abut  one  another  as  a  singe  refinement 
region  sounded  by  a  single  ghost  cell  zone,  2)  uses  ghost  cells  strictly  to  provide  boundary  data  to 
the  refinement  region,  3)  uses  Hermite  interpolation  to  construct  third  order  in  time  ghost  cell  data 
for  intermediate  time  levels  which  arise  in  order  to  satisfy  the  CFL  of  the  fine  mesh  in  the  refine¬ 
ment  region,  4)  uses  WENO  interpolation  in  space  to  provide  an  initial  a  5th  order  approximation 
to  the  course  soliton  on  the  fine  mesh  and  to  update  the  ghost  cell  regions,  5)  uses  the  WENO  indi¬ 
cator  to  disced  were  to  refine,  6)  the  method  use  a  3rd  order  TVD  RK  method  to  do  time  integration. 

We  have  show  that  the  method  is  indeed  5th  order  in  space  and  3rd  order  in  time  for  smooth  prob¬ 
lems.  Eurther,  we  have  shown  5rd  order  AMR  is  far  less  diffusive  than  low  order  AMR,  requiring 
less  refinement  because  it  is  less  diffusive,  for  a  range  of  test  problems.  This  include  the  ID  dam 
break  problem,  the  ID  blast  wave  problem,  and  the  2D  double  Macoh  refection  problem.  In  all 
cases,  the  method  give  superior  results,  coasting  less  computational  effort  for  a  given  accuracy  than 
the  low  order  AMR  method. 

We  are  working  to  extend  these  ideas  to  a  constrained  transport  solver  for  MHD  and  to  extend 
these  ideas  to  a  high  order  Vlasov  solver. 
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